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Abstract 


The coupled exciton-vibrational dynamics of a 3-site Fenna-Matthews-Olson (FMO) model 
is investigated using the numerically exact multilayer multiconfiguration time-dependent Flartree 
approach. Thereby the vibrational mode specific coupling to local electronic transitions is 
adapted from a discretized experimental spectral density. The solution of the resulting time- 
dependent Schrodinger equation including three electronic and 450 vibrational degrees of free¬ 
dom is analyzed in terms of excitonic populations and coherences. Emphasis is put onto the 
role of specific ranges of vibrational frequencies. It is observed that modes between 160 and 
300 cm~ 1 are responsible for the subpicosecond population and coherence decay. Further, it 
is found that a mean-field approach with respect to the vibrational degrees of freedom is not 
applicable. 


Introduction 

The discovery of long-lasting coherent oscillations in two-dimensional spectra obtained for a num¬ 
ber of photosynthetic antenna complexes up to physiological temperatures 13 has triggered a host 
of theoretical investigations (for a review, see ref 0). A controversial discussion focussed on the 
role of vibrational degrees of freedom (DOFs). Traditionally, vibrations have been considered to 
form a heat bath, thus leading to phase and energy relaxation and serving the purpose of directed 
downhill energy transfer. 5 6 With the decoherence time of the observed oscillations considerably 
exceeding the expectations based on typical line widths of electronic transitions, the coupling to 
vibrations was reconsidered with the result that it might be more important than anticipated based 
on the smallness of the Huang-Rhys (HR) factor. 7 8 In fact a special role was ascribed to a vi¬ 
brational mode having a frequency that matches the electronic energy gap between neighbouring 
bacteriochlorophyll sitesJ 7 9 10 Such a mode, which for the relevant sites one to three in the FMO 
complex should have a frequency of about 180 cm -1 , could facilitate vibrationally-assisted reso¬ 
nance transfer. In terms of the eigenstates of the Coulomb-coupled excitonic and vibronic excita¬ 
tions there might be a considerable quantum state mixing 1112 making an analysis of the observed 
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spectra in terms of pure electronic or vibronic excitations difficult. 13 14 Interestingly, a popular 
experimentally determined spectral density features a peak in the 180 cm -1 region. 15 

Simulations of coupled exciton-vibrational dynamics are usually performed using non-pertur- 
bative and non-Markovian methods such as the hierarchy equation of motion 16 18 or path integral 
approachesJ- 19 Here, the dynamics is obtained explicitly only for the relevant system, whereas 
the bath is traced out in the reduced density operator. 20 Apart from a few examples, where a 
single vibrational mode has been kept alive in the relevant system, 13 21 22 the latter consists of the 
excitonic DOFs only. Hence, there is no access to the explicit dynamics of the vibrational modes 
within these system-bath approaches. In passing we note that the inclusion of a single vibrational 
mode into the relevant system has, of course, some arbitrariness with respect to the choise of the 
modes’ properties. 

In the present contribution we study the coupled exciton-vibrational dynamics from a com¬ 
pletely different perspective. Specifically, we solve the time-dependent Schrodinger equation for 
a high-dimensional model of the FMO complex. This allows us, for the first time, to access infor¬ 
mation on the dynamics of electronic ground/excited state vibrational/vibronic excitations during 
excitation energy transfer in this complex. This becomes possible due to the special structure of 
the Frenkel exciton Hamiltonian, which is ideally suited for the combination with the multilayer 
multiconfiguration time-dependent Hartree (ML-MCTDH) wave packet propagation method. 23 26 


Theory 

We will use the standard Frenkel exciton Hamiltonian to describe an aggregate with N agg sites 
(site index m), each site having the excitation energy E m , and different sites being coupled by the 
Coulomb interaction J mn : 

^agg 

Hex = (&mnEm T Jmn ) \ m ) (^| • (1) 

m,n=l 

Local electronic states are restricted to the ground, | g m ), and excited, \e m ), states, i.e. the local 
or diabatic one-exciton states are given by | m) — \e m ) YJ n ^ m I gn)- For the present simulations we 


3 






will adopt the 8-site FMO Hamiltonian reported by Moix and coworkers. 27 It is reduced to a 3-site 
model (sites 1-3) as suggested by the simulations of population dynamics in refl27l(for full set of 
parameters see Supplementary Information (SI)). Indeed, starting from an initial population of site 
one only, the populations of individual sites other than 1-3 never exceed about 12% within the first 
1 ps. This finding is supported by the decomposition of eigenstates into the local basis for full and 
reduced models as given in Table SI (SI). In the 3-site model the exciton Hamiltonian matrix is 
given by (in cm -1 ) 


H„ = 


( 310 -98 6 ^ 


-98 230 30 
6 30 0 


( 2 ) 


Diagonalization of this matrix yields the one-exciton eigenstates, \a), whose decomposition 
into the local states | m) is given by a) = C a . m \ m ) and shown together with the eigenenergies 
in Figure |2j The local vibrations at site m are described in harmonic approximation by the set 
of dimensionless normal mode coordinates {Q m ^} with frequencies s}, i.e. the vibrational 
Hamiltonian reads 


ffv«, = E£ 

m 


fr(Q, 


'm,% 




(3) 


Here, the notation E, G m indicates the summation over all vibrational modes of site m. Exciton- 
vibrational coupling is accounted for within the linearly shifted oscillator (Huang-Rhys) model 
given by 

^ex-vib = II ft (Omt j2S m £ Q m £ I m) (m\ . (4) 

tn 


The coupling of a particular mode to the local electronic transition and thus the degree of vibronic 
excitation, is characterized by the HR factor S m £. In system-bath models this coupling is typically 
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expressed via the spectral density^ 20 


J m {(o)=A £ S m ^8(co-co m ^), (5) 

i,6m 

where A is a constant to be specified below. Note that J m ((o) corresponds to the spectral density 
for monomeric bacteriochlorophyll a in the actual FMO protein environment. There are several 
calculations of this spectral density based on classical molecular dynamics and taking into account 
the protein and solvent environment. 28 31 However, the reported results differ considerably, with 
respect to their structure and integrated intensity; the latter gives the total HR factor for site m: 
Stot = A 1 fd(oJ m (co) — L^nSm.^ 20 Therefore, we will use the experimentally determined spec¬ 
tral density of Wendling et al. 15 shown in Figure [ 2 ] It has been obtained from low-temperature 
site-selected fluorescence, measured for the energetically lowest pigment of the complex under the 
assumption that the Coulomb coupling to the other pigments can be neglected. The total HR factor 
was determined as 0.42. Note, however, that it has recently been pointed out that the correction to 
the HR factor of particular modes due to the Coulomb coupling can be substantial, i.e. for modes 
in the 180 cm -1 region up to a factor of 1.5. 12 Notice further that the experimental spectral density 
covers the range up to about 350 cm -1 only, whereas some calculations predict distinct structure 
up to 2000 cm -1 . 28 However, excitonic transition energies between strongly coupled pigments 
are essentially located in the range up to 300 cm -1 , what justifies to neglect higher frequency 
modes. In rcfUTla detailed comparison of this low-frequency part of 7 (cd) with results of Quantum 
Mechanics/Molecular Mechanics calculations was given. Although different in details, the calcu¬ 
lations could in particular be used to justify the assumption of the HR model leading to eq ?? and 
eq ??. In the present model we will use the original spectral density from ref [15] and discretize it 
into 150 modes within the interval [2 : 300] cm -1 . The amplitudes of the individual HR factors 
have been adjusted homogeneously via the constant A in eq ?? such as to preserve 5 to t = 0.42 upon 
summation. 

The time-dependent Schrodinger equation will be solved employing the ML-MCTDH method 
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(for a review, see refl32l). Here, the the state vector is expanded in terms of the local exciton basis 
according to 


l^(Q;0>=I>m(Q;0l™> ) (6) 

m 

where the nuclear coordinates are comprised into the D = N agg x = 3 x 150 = 450 dimensional 
vector Q. The nuclear wave function for each diabatic state is expanded into MCTDH form 

n i\- n i D 

Xa(Q,t)= £ c'“>(7) 

jl-jD 

(ex') 

Here, the C^ ; U) (l ) are the time-dependent expansion coefficients weighting the contributions 
of the different Hartree products, which are composed of nj k single particle functions (SPFs), 
< P^\Qk;t ), for the Ath degree of freedom in state a. In ML-MCTDH the SPFs itself describe 
multi-dimensional coordinates that are expanded into MCTDH form. 25 2633 This nested set of 
expansions can be represented by so-called ML-MCTDH trees. 26 The particular choice of this tree 
may have a strong influence on the required numerical effort. 4 In the following simulations we use 
a grouping according to the magnitude of the HR factor as detailed in the SI. There we also give 
the other parameters of the ML-MCTDH setup. Wave packet propagations have been performed 
using the Heidelberg program package. 34 Temperature effects due to the thermal population of 
vibrational states in the electronic ground state are not included. In all case the initial conditions has 
been a vertical Franck-Condon transition at site 1 and the propagation time was 1 ps. Convergence 
of the ML-MCTDH setup has been monitored by means of the natural orbital populations] 24 In the 
full/reduced model the largest population of the least occupied natural orbital was ~ 10 -4 /~ 10 6 . 

Results and Discussion 

Figure [T] shows the diabatic population and coherence dynamics of the 3-site model. The diabatic 
populations show a clear beating between sites 1 and 2, overlaid by a damping of both oscillation 
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Figure 1: Upper panel: Diabatic population dynamics for the 3-site FMO model. Lower panels: 
Coherence dynamics as indicated in terms of the off-diagonal elements of the electronic density 
matrix in the diabatic, i.e. site-based, representation. Note that the long-time behaviour of the 
coherences reflects the mixing of the adiabatic states with respect to the diabatic ones (cf. Figure [2j). 


amplitudes. For site 3 one observes only rather small amplitude oscillations on top of an other¬ 
wise monotonously increasing population. A Fourier decomposition of the population dynamics 
yields a dominant component at 220 cm -1 (period 151 fs). This does not agree with the bare elec¬ 
tronic oscillation period, which would be 160 fs (208 cm -1 ). It is this difference which actually 
reflects the modified energy level structure of the coupled exciton-vibrational system. In principle 
this oscillatory behaviour is qualitatively similar to previously published results. 27 35 36 There are 
quantitative differences due to different Hamiltonian parameters 35 36 and spectral densities. 27 35 36 
In addition, the present model assumes a discretized bath and does not include temperature, which 
might be the reason for the population oscillations to last over the whole time interval. 

The coherences, p m „(t ) = (m| v F(t)) (*F(/)| n), between the different sites also show damped 
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oscillations as seen in the lower panels of Figure |T] The oscillations of coherence density matrix 
elements have been related to the off-diagonal peaks in the two-dimensional spectroscopy of the 
FMO complex. 2 37 Note, however, that the actual assignment to specific pairs of exciton eigenstates 
depends on the parametrisation of the Hamiltonian. Experimentally, in ref[2] oscillations at about 
160 cm -1 have been reported, irrespective of temperature in the range between 77 K and 150 K. 
In the present case, which formally corresponds to zero temperature, the dominant frequencies in 
this range are at 220 cm -1 (pn) and at 190 cm -1 (P 23 ); both frequencies occurring for p 13 . We 
also observe that the amplitudes of these oscillations decay considerably within the given time 
window due to the Coulomb coupling induced mixing between electronic and vibronic/vibrational 
excitations. Apparently, the finite but large number of vibrational DOFs is sufficient to mimic 
excitonic energy relaxation and coherence dephasing. 



I I I I I ~~I- 1^0 

0 0.5 1 1.5 2 2.5 3 diab. state 


Figure 2: Vibronic (upper row) and vibrational (lower row) excitation dynamics of the 3-site FMO 
model using a discretization of the spectral density into 150 equally spaced modes in the interval 
[2 : 300] cm -1 for each monomer. Note that the vibrational excitation of site 3 is close to zero 
and not shown. The spectral density is given in the left panels. The lower right panel shows the 
spectrum of one-exciton eigenstates and their decomposition into amplitudes of (local) diabatic 
states. The color code covers the interval from 0 to 3 cm -1 ; zero point energy has been subtracted. 
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In order to investigate the vibronic/vibrational excitation more closely we consider the dynam¬ 


ics in the different potential energy surfaces. Specifically, we define the operator 



whose expectation value gives the vibrational energy in the electronic ground states for monomer 
m, no matter the electronic population of the other monomers n ^ m. 

The vibronic energy at monomer m follows from the expectation value of 



Figure [2] shows the dynamics of vibronic and vibrational excitations for the situation of Fig¬ 
ure [I] First, we notice that the degree of excitation is rather small. Expectation values of the energy 
do not exceed 3 cm -1 for individual modes. Second, there is a marked difference between the in¬ 
dividual sites and the dynamics in their respective potential energy surfaces. This can be explained 
either in a dynamical or an eigenstate picture. Within a dynamical picture vibrational excitation in 
the electronic ground state corresponding, e.g., to the site of initial electronic excitation, is a conse¬ 
quence of the finite transfer time between sites 1 and 2. If the latter is longer than the time scale of 
vibrational motion the initial vibronic wave packet moves out of the original Franck-Condon win¬ 
dow. Thus, upon deexcitation during the transfer event, vibrationally excited states are populated 
in the electronic ground state. We emphasise that the oscillations due to exciton transfer observed 
in Figure |T| are reflected as well in the vibrational and vibronic excitations. Eventually, the exciton 
population is accumulated at site 3 and so is the vibronic excitation. The trapping at site 3 comes 
along with negligible vibrational excitation at this site (not shown). However, the nonequilibrium 
vibrational excitation at sites 1 and 2 is maintained even at 1 ps. Note that in order to describe vi¬ 
brational relaxation, which would be relevant at longer time scales, one should include anharmonic 
interactions among the harmonic oscillator modes. 

Third, considerable excitation is restricted to certain mode frequency ranges and here most 
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notably between about 160 to 200 cm -1 . In the dynamical picture this could be explained by 
the fact that the electronic time scale is about 160 fs, i.e. only modes with higher frequencies 
have a chance to leave the Franck-Condon window and, therefore, become vibrationally excited. 
The selectivity of strong vibronic excitation in particular at site 3 is better explained in terms 
of a resonance effect, i.e. vibrationally-assisted exciton transfer. This is similar to the vibronic 
enhancement of exciton transport discussed for allophycocyanine dimers in refl38l In Figure [2] we 
also show the electronic eigenenergies and the decomposition of eigenstates into the local exciton 
basis. The energetic separation between adjacent eigenstates is around 200 cm -1 , which explains 
the preference of excitation for this region of the spectral density. 



Figure 3: Vibronic (upper row) and vibrational (lower row) excitation dynamics of the 3-site FMO 
model using a reduced description where only vibrational modes with frequencies in the interval 
[159 : 300] cm -1 were included. The lower right panel compares the diabatic populations for full 
(dots) and reduced (lines) models. The color code covers the interval from 0 to 3 cm -1 ; zero point 
energy has been subtracted. 


In order to investigate the role of the high frequency part of the spectral density closer, we 
have studied a reduced model where only modes in the range [159 : 300] cm -1 were included. The 


10 
















resulting dynamics is shown in Figure [3j The lower right panel compares full and reduced models 
in terms of the diabatic population dynamics. Clearly, the differences between the two models are 
small and can be found mostly in the amplitudes of the oscillations, which are generally larger in 
the reduced model. 

The vibrational and vibronic excitation also behaves similar, with the main difference being 
higher excitation energies obtained for the reduced model which reaches 4 cm -1 per mode. Com¬ 
paring Figure [2] with Figure [3] this can be interpreted as the effect of vibrational energy redistribu¬ 
tion, which is more pronounced in the full model due to the larger spectral range that is available. 
We emphasise again that for a single site the vibrational DOFs are decoupled by construction of the 
model. Vibrational energy redistribution is possible solely due to the Coulomb coupling between 
the sites. 



Figure 4: Vibronic (upper row) and vibrational (lower row) excitation dynamics of the 3-site 
FMO model using the TDH approximation, eq ??. The lower right panel compares the diabatic 
populations for the ML-MCTDH (dots) and TDFI (lines) models. In contrast to Figure [2] and 
Figure [ 3 ] the color code covers the interval from 0 to 1.5 cm -1 only; zero point energy has been 
subtracted. 
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The ML-MCTDH approach provides a many-body wave function which is exact in the sense 
of the given numerical convergence thresholds. One might ask to what extent correlations actually 
matter. After all the isolated monomer is described by uncoupled vibronic excitations within the 
present model. The effect of Coulomb coupling on the correlations between vibrational DOFs is 
most straightforwardly accessed by comparison with the time-dependent Hartree (TDH) approxi¬ 
mation to the wave function in eq ?? which reads 

*« TDH) (Q, t) = C (a) O)0i (a) (Gi; 0 • • • <t>D ] (Qd; 0 • (10) 

In passing we note that the CPU time is reduced by a factor of about 600 when using eq ?? 
instead of eq ??. However, the resulting dynamics shown in Figure [4] is not at all close to the 
ML-MCTDH simulation in Figure [2] We start by comparing the population dynamics. Not only 
that the oscillation period for population exchange between sites 1 and 2 changes, the exciton 
population also stays at sites 1 and 2 for a much longer time in TDH, i.e. the population trapping 
at site 3 is not very efficient within the first 1 ps. This can be attributed to lack of sufficient 
vibrational energy redistribution, which comes along with the less flexible TDH wave function. 
This is further supported by the vibrational and vibronic energy distributions shown in Figure |4j 
First, the degree of excitation is much smaller than in Figure [2} notice the different scales in the 
two figures. Second, apart from the oscillation pattern which is imprinted by the long-lasting 
oscillations of the populations, the spectral distribution of vibrational and vibronic excitations is 
also modified. Thus we conclude that the TDH approximation is not suitable for the description of 
the present FMO model. 

Conclusions 

In summary, we have presented a numerically exact quantum dynamics simulation of the Coulomb- 
coupled exciton-vibrational dynamics of a 3-site model of the FMO complex at zero temperature. 
This has become possible by combining the simplicity of the Frenkel exciton/linearly displaced 
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harmonic oscillator Hamiltonian with the ML-MCTDH approach, thereby incorporating vibronic 
coupling via an experimentally determined spectral density. Adopting this model, numerically ex¬ 
act refers to the used convergence threshold for the natural orbital populations within ML-MCTDH. 

Most notably we have found that typical time scales for excitonic energy and phase relaxation 
can be reproduced by a finite discretization of an experimental spectral density. Further, in accor¬ 
dance with previous suggestions 7 39 modes in the 180 cm -1 region play an important role for the 
transfer. The novel feature of explicit full-dimensional quantum dynamics used in the present sim¬ 
ulation, however, is that it gives access to the dynamics of all DOFs. This allowed to discover that 
a large number of modes between 160 and 300 cm -1 are actually excited in the electronic ground 
and excited states. This finding points to a possible dimensionality reduction in future studies, 
which may facilitate to account for all eight sites of the FMO. In other words, the reduction to a 
single mode represents a rather strong approximation. However, a simplification in terms of the 
mean-field description as provided by the TDH approximation is not suitable as it gives qualita¬ 
tively different results. In other words, there is a strong correlation in the many-body dynamics of 
excitonic and vibrational DOFs. We conclude by pointing out that the present simulation provides 
further evidence for the exceptional role that is played by vibrational and vibronic excitations in 
photo synthetic energy transfer. 40 
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High-dimensional quantum dynamical simulations are performed for a FMO model which is based 
on an experimental spectral density (left panel). In the right panel the vibronic excitation is shown 
for the terminal site of the complex. 
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